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OO ■ A new variational principle for optimizing thermal density matrices is intro- 

^_^ ' duced. As a first application, the variational many body density matrix is written 

fH , as a determinant of one body density matrices, which are approximated by Gaus- 

O , ' sians with the mean, width and amplitude as variational parameters. The method 

I ^ is illustrated for the particle in an external field problem, the hydrogen molecule 

S . and dense hydrogen where the molecular, the dissociated and the plasma regime 

C/3 ' are described. Structural and thermodynamic properties (energy, equation of state 

^ , and shock Hugoniot) are presented. 

c/j ■ 

ZJ , I. INTRODUCTION 

C/2 . 

^^' Considerable effort has been devoted to systems where finite temperature ions (treated either 

1-^ I classically or quantum mechanically by path integral methods) are coupled to degenerate electrons 

^_j- on the Born-Oppenheimer surface. In contrast, the theory for similar systems with non-degenerate 

electrons (T a significant fraction oiTpermi) is relatively underdeveloped except at the extreme high 

^-H , T limit where Thomas-Fermi and similar theories apply. In this paper we present a computational 

^ ' approach for systems with non-degenerate electrons analogous to the methods used for ground 

^Q , state many body computations. 

^^ ' Although an oversimplification, we may usefully view the ground state computations as 

^—v I consisting of three levels of increasing accuracy M. At the first level, the ground state wave 

^_^ . function consists of determinants, for both spin species, of single particle orbitals often taken from 

^\ ' local density functional theory 

o>: 

o : *Gs(R) 

>>: 

,S^ , The majority of ground state condensed matter calculations stop at this level. 

Mh' If desired, additional correlations may be included by multiplying the above wave function by 

^ I a Jastrow factor, Yii /(^ii)i where the / will also depend on the type of pair (electron-electron, 

l^ • electron- ion) . Computing expectations exactly (within statistical uncertainty), with this type of 

r> I wave function now requires Monte Carlo methods. 

j^ ' Finally diffusion Monte Carlo ^,^ methods using the nodes of this wave function to avoid 

■ ■ ■ ' the Fermion problem may be used to calculate the exact correlations consistent with the nodal 

structure. 

The finite temperature theory proceeds similarly. Rather than the ground state wave function 
a thermal density matrix 

p(R, R'; /3) = (R|e-'3«|R') = ^ e-'5^=M/,(R)vI/,(R') (2) 

s 

is needed to compute the thermal averages of operators 



$i(ri) ... $jv(ri) 

$l(rAr) ... $Ar(rAr) 



(1) 



At the first level, this many body density matrix may be approximated by determinants of 
one-body density matrices, for both spin types, as well as the ions 



p(R,R';/3) 



pi{ri,r[;p) ... pi{rN,r[;P) 
Pi{ri,r'f^;/3) ... pi{rN,r'^;P) 



(4) 



The Jastrow factor can be extended to finite temperature and the above density matrix 
multiplied by Yii ,- f{fiji''''ipP)- In particular, the high temperature density matrix used in path 
integral computations has this form. 

Finally, the nodal structure from this variational density matrix (VDM) may be used in 
restricted path integral Monte Carlo (RPIMC) p|-p|. This method has been extensively applied 
using the free particle nodes. One aim of the present work is to provide more realistic nodal 
structures as input to RPIMC. 

This paper considers the first level in this approach. The next section is devoted to a general 
variational principle which will be used to determine the many body density matrix. The principle 
is then applied to the problem of a single particle in an external potential and compared to exact 
results for the hydrogen atom density matrix. After a discussion of some general properties, 
many body applications are considered starting with a hydrogen molecule and then proceeding to 
warm, dense hydrogen. It is shown that the method and the ansatz considered can describe dense 
hydrogen in the molecular, the dissociated and the plasma regime. Structural and thermodynamic 
properties for this system over a range of temperatures (T= 5 000 to 250 000 K) and densities 
(electron sphere radius Vg = 1.75 to 4.0) are presented. 



II. VARIATIONAL PRINCIPLE FOR THE MANY BODY DENSITY MATRIX 

The Gibbs-Delbruck variational principle for the free energy based on a trial density matrix 

F < TrlpH] + kT Tr[/5 In p] (5) 

where 

p = p/Tr[p] (6) 

is well known and convenient for discrete systems (e.g. Hubbard models) but the logarithmic 
entropy term makes it difficult to apply to continuous systems. Here, we propose a simpler 
variational principle patterned after the Dirac-Frenkel-McLachlan variational principle used in 
the time dependent quantum problem p]. Consider the quantity 

as a functional of 

I{Q) = Tr{Q + npf (9) 

with p fixed. / (O) = when Q satisfies the Bloch equation, 6 = —Tip, and is otherwise positive. 
Varying / with gives the minimum condition 



Tr [Se{e + Hp)] =0 . (10) 

This may be written in a real space basis as 

f I (5e(R', R; (3) [e(R, R'; (i) + Hp(R, R'; (3)] dRdR' = (11) 

or, using the symmetry of the density matrix in R and R', 

I I (5e(R, R'; /3) [e(R, R'; (3) + H/9(R, R'; (3)] dRdR' = . (12) 

Finally, we may consider a variation at some arbitrary, fixed R' to get 

I (5e(R, R'; (3) [e(R, R'; (3) + np{R, R'; /?)] dR = VR'. (13) 

It should be noted that in going from Eq. O to Eq. O a density matrix symmetric in R and 
R' is assumed, which is a property of the exact density matrix. If the variational ansatz does not 
manifestly have this invariance Eq. n3 minimizes the quantity, 

/"[e(R,R';/3)+7^p(R,R';/3)]^dR = . (14) 

We propose solving this equation by parameterizing the density matrix with a set of parameters 
Qi depending on imaginary time f3 and R', 

p(R, R'; /?) = p(R, gi, . . . , g™) where g,(R'; /3) (15) 



so 



e(R,R,/3) = L^ 9^ = ^* % • ^''^ 

I— 1 z— 1 

In the imaginary time derivative B only variations in q and not q are considered since p is fixed so, 

5e(R,R';/?)=f;<5g,(R';/3)^4r^ ' ^^^^ 

Using this in equation 03 gives for each variational parameter, since these are independent. 



/-^(e + Hp)dR = o 

J dqj 



(18) 



This reveals the imaginary-time equivalent to the approach of Singer and Smith pO| ] for 
an approximate solution of the time dependent Schodinger equation using wave packets (see 
section 



[II). Introducing the notation 



.. - ^ (19) 

and using Eq. Hq, the fundamental set of first order differential equations for the dynamics of the 
variation parameters in imaginary time follows from Eq.. Hq as, 

Pj pHp dR + ^ q, jpj p, p^dR = (20) 



or in matrix form 



where 



and the norm matrix 



with 



IdH ^ • 

^^+AA.^^O (21) 



M, 



H = / pTip rfR (22) 



yp,,,p^.R=lhn^^-^ (23) 



N = I p(R, g ; /3) p(R, q' ; f3) dR . (24) 

The initial conditions foUow from the free particle limit of the density matrix at high temperature, 

p{R, R'; p) -> exp [-(R - R') V4A/3] /{AnXpf^/^ where A = h^ /2m . (25) 

Various ansatz forms for p may now be used with this approach. After considering the analogy to 
real time wave packet molecular dynamics, the principle is first applied to the problem of a particle 
in an external field. 



III. ANALOGY TO REAL-TIME WAVE PACKET MOLECULAR DYNAMICS 

Wave packet molecular dynamics (WPMD) was first used by Heller mU and later applied to 
scattering processes in nuclear physics [|l2[ and plasma physics [ [l3|jl^ ]. An ansatz for the wave 
function ijj{qu) is made and the equation of motions for the parameters Qi, in real time can be 
derived from the principle of stationary action |12[| , 

sfdtL^O , L{q,{t),q,{t))^{'4j\idt~rL\^) (26) 

This leads to a set of first order equations, which provides an approximate solution of the 
Schrodinger equation. However, this principle cannot be directly applied to the Bloch equation 
because there is no imaginary part in the density matrix. For this reason, we followed in our 
derivation in section O the principle of Dirac, Frenkel and McLachlan H, which minimizes the 
quantity 

f \n^p - ihe\^ dt, ^ = ^- (27) 

This method was employed in pJ]| to obtain the dynamical equations in real time. 

The VDM approach and WPMD method share the zero temperate limit, which is given by 



the Rayleigh-Ritz principle (see section VA). At high temperature, the width of wave packets 
in WPMD grows without limits, which is a known problem of this method ||l^,|l6[. In the VDM 
approach, the correct high temperature limit of free particles is included. The average width shown 
in Fig. |lO| can be used to verify the attempts to correct the dynamics of the real time wave packets 



IV. EXAMPLE: PARTICLE IN AN EXTERNAL FIELD 



As a first example, we apply this method to the problem of one particle in an external potential 

n = -XV^ + V{r) . (28) 

The one-particle density matrix will be approximated as a Gaussian with the mean m, width w 
and amplitude factor Z?, 



pi(r,r',/3) = (7r«;)-3/2 



exp 



1 



{r-mf +D 



(29) 



as variational parameters. The initial conditions at /3 — > are w — 4A/3, m = r' and D = in 
order to regain the correct free particle limit, Eq. G^ For this ansatz H, defined in Eq. E3 as 



H^ fpHpdr= (^ + ^f°l 



o2I> 



where 



and 



(27rw)3/2 

yln] ^ ( J_)3/2 /"(j. „ m)"l^(r)e-'('^-™)'/"'dr 
nw J 



N= pp'dr= [7r(w + ?«')]"^/2gxp{-(m-m')V(«^ + if'')}exp(D + i:'') 
From Eq. ^, the equations for the variational parameters are. 



m 

b 



-2V[il 



Vm _ ll/[2] 



(30) 



(31) 



(32) 



(33) 
(34) 
(35) 



In absence of a potential, the exact free particle density matrix is recovered. The harmonic oscillator 
case is also correct since the Gaussian approximation is exact there. For a hydrogen atom, A = 1/2, 
V{r) = — 1/r and 



F™ = erf |m 

m 



y[l] ^ iii^ 

-m? 4 



v/27^) 
' I my2/w ] 



erf ( my2lw — \ e 

^ / V TTW 



-2m^/i 



y[2] ^ /^e-2" /» + — yM 
V 27r 4 



(36) 
(37) 

(38) 



At low temperature, the density matrix as a function of r goes to the ground state wave function as 
discussed in more detail in the next section. One expects this to be a fixed point of the dynamics 
of the parameters m and w determined by in = and w — while D — —Eq. The (3 ^ 00 fixed 
point: m = Q, w ^ 9tt/8, D = 4/37r (atomic units) corresponds to the well known Rayleigh-Ritz 
variational result for a Gaussian trial wave function 



3/2 
*o(r) = ( ^ ) exp(~8rV9^) • 



(39) 



In ground state variational studies, addition of two more Gaussians brings the ground state energy 
to within 0.6% of exact and similar improvement would be obtained here. 

Results at finite (3 require a numerical solution, which is illustrated in the figure below 
comparing the Gaussian variational density matrix with the exact p^ and the free particle density 
matrix at several temperatures for the initial condition r' = 1. At high temperatures (/? = .05 and 
(3 = .25) the Gaussian approximation correctly reproduces the limiting free particle density matrix. 
At lower temperatures, the cusp in the exact density matrix due to the Coulombic singularity at 
the proton becomes evident and the peak shifts to the origin somewhat faster than the Gaussian 
variational approximation. As P increases the exact result grows faster than the variational since 
the correct energy, -0.5, is lower than — 4/37r but the Gaussian variational approximation remains 
rather accurate for r > 1. The free particle density matrix remains centered at r = 1 and beyond 
(3 — 0.5 (T = 54.4 eV) bears little resemblance to the correct result. 




FIG. 1. Comparison of the Gaussian variational approximation (circles) with the exact density matrix 
p(r, r';/3) (solid line) for a hydrogen atom. The free particle density matrix (dashed line) is also shown. 
The plotted r is along the line from the proton at the origin (marked by the vertical bar) through the 
initial electron position r' = 1. 



V. VARIATIONAL DENSITY MATRIX PROPERTIES 



A. Zero Temperature Limit 



In the preceding section, it was shown that for the hydrogen atom the Gaussian variational 
density matrix, as a function of R converges at low temperature to the Gaussian ground state wave 
function given by the Rayleigh-Ritz variational principle. It is generally true that the Rayleigh-Ritz 
ground state corresponds to a /3 — > oo of the variational density matrix as we now show. 



The Rayleigh-Ritz principle states that for any real parameterized wave function 
^(R, qi, . . . ,qm) the variational energy 

_ /V;(R)H^(R)dR 

is greater than or equal to the true ground state energy even at the minimum determined by 

±-^E{{q})^Qyk. (41) 

For the VDM ansatz, an amplitude parameter D is assumed such that 

p(R, R'; P) = e^^^'^/^^plR, {g(R'; /?)}) . (42) 

As in the one particle example, it is expected that at low temperature, /3 — > oo, the other q^ — + 
while D -^ constant. From this assumption, Eq. |2^ implies that as /? ^ oo 

f.i.f.0 (43) 

oqk oqk 

for all variational parameters, where wc have defined H = J pTip dH and N = J p^ rfR. Since 
dH/dD = 2H and dN/dD = 27V, Eq. || for qt = D imphes D = -H/N = -^o so Eq. || may 
be rewritten as 

at the /? ^ oo fixed point. With the correspondence 

p(R,{g(R',/3)})^e^(^'^'5)^(R,M) , (45) 

this is equivalent to Eq. |4| and thus the Rayleigh-Ritz ground state corresponds to a zero 
temperature fixed point in the dynamics of the parameters. 

Z? is a function of R' and /3, which is calculated by integrating from /3 = with Eq. g5| as 
initial conditions. The zero temperature limit of Z? is a constant, —Eq, which means in the low 
temperature limit D can written as 

D(R';/3) = -(3Eo + .f(R') . (46) 

The function /(R') can be rewritten as, 

/(R') = ln{^o(R')[l + '5(R')]} , (47) 

where the function (5(R') is introduced to describe the variational error in the solution of the Bloch 
equation. It is identical to zero if the variational ansatz includes the exact solution. It leads to 
loss of symmetry in R and R', which will discussed in the next section. Eq. ^ now reads, 

p(R, R', /? ^ (») - e-'5^°^o(R)^o(R') [1 + Sin')] (48) 

For certain potentials, several fixed points of the dynamics can exist. From Eq. ^, it follows 
that only the lowest energy state contributes to physical observables calculated from Eq. |[ This 
completes the argument that the zero temperature limit of the VDM correspond to the Rayleigh- 
Ritz ground state. 

In case of an anti-symmetrized ansatz for the density matrix, one can show that the fixed 
point of the dynamics in imaginary time corresponds to the Rayleigh-Ritz ground state for an 
anti-symmetrized wave function. 



B. Loss of Symmetry 

The exact density matrix is symmetric under R <-> R'. Since we have singled out R' as the 
initial point for the imaginary time dynamics, it is not clear that the approximation given in Eq. E9^ 
automatically satisfies this condition. For the free particle limit and the harmonic oscillator, where 
the Gaussian is the exact solution, it obviously does but in general it does not. 

As a specific example, consider again the ground state limit of the hydrogen atom where the 
Gaussian VDM approximation. Eq. E9 then reads, 



lim p(r,r';/3) = e^(''''« (8/97r^) 



2>|3/2 -8r^/97r 



(49) 



For this to be symmetric under r <-> r', we must have 

lim D{r';P) = -8/ VStt + c(/3) 



(50) 



and from the result for D, lim^j^oo c{f3) = 4(3 /Sir + c\. 

Figure compares the D{r, 0) from the Gaussian VDM with Eq. 

c(/3) =4/3/37r + 3/21n2. 
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FIG. 2. D{r,P) from the Gaussian approximation in the ground state limit (solid line) of the hydrogen 
atom. Deviations of this function from linearity indicate a breakdown of symmetry in the Gaussian 

„2, 



approximation for o(r, r 
ground state Eq. pS. 



; /3) . The dashed line is -8r /97r+4/3/37r+3/2 In 2 expected from the Rayleigh-Ritz 



There are several consequences of this small violation of R <~> R' symmetry. As shown generally 
in the section above, in the /3 ^ oo limit —D is the Rayleigh-Ritz variational ground state energy 
for a Gaussian wave function, which for the hydrogen atom is Eq = — 4/37r = —0.4244. Because of 
the loss of symmetry this is not the same as the energy given by the estimator 



(E) = (H) 



Tr[7^p] 
Tr[p] 



(51) 



in the /3 — > oo limit, which for the hydrogen atom gives the more accurate result (E) = —0.4709. 
This will be seen again below for the hydrogen molecule where Eq. M also gives more accurate 



ground state energies. Other consequences are less pleasant. Although the energy is more accurate 
the virial theorem, {K) — — (U) / 2, between the kinetic and potential energy is violated by about 
3% (while both are more accurate than the usual ground state variational Gaussian result). This 
has consequences for calculating the equation of state particularly at low density. Slightly more 
complicated, explicitly symmetric forms for the VDM could be used but in this paper we will 
continue to explore the basic Gaussian approximation. 



C. Thermodynamic Estimators 

Since the VDM, except in the simplest cases, is not exact various estimators for the same 
quantity will differ. For example the variational principle introduced in section II consists 
essentially in globally minimizing the squared difference between dp/d(3 and Tip, either of which 
can be used in estimating the energy. As mentioned above the energy estimator Eq. |5l| and its 
kinetic and potential energy pieces do not automatically satisfy the virial theorem for Coulomb 
systems at low density. As an alternative to Eq. ^, one can use the thermodynamic estimators, 

W=-(|jl»p). (52) 

{r)^-|(A,„,). (53) 

for the total, kinetic and potential energy. These estimators satisfy 

(E) = (T> + {V) (55) 

by the following argument. Any function / = /(/3A,/3e^) satisfies 

From Eq. ^it follows that all parameters qi — qi(R.'; P, A, e^) have this property and therefore so 
does the variational density matrix. 

In the zero temperature limit, the thermodynamic estimators satisfy the virial theorem, which 
is also satisfied by any exact and any variational Rayleigh-Ritz ground state. From the zero 



temperature limit of the VDM given by Eq. Wq and the 1//3 factor in Eqs. 53 and p4, it is seen 
that the symmetry error (5(R') is unimportant in this limit. It should be noted that calculating 
the derivatives for (T) and {V) increases the numerical work. The pressure is estimated from 

3 {P)v^2 (K) + {V) . (57) 



VI. MANY PARTICLE DENSITY MATRIX 

We represent the many particle density matrix by a determinant of one-particle density matrices 
(Eq. |4|). It can written as, 

p(R,R',/?) = ^epn'°i('^^-'^k'/5) = E^''^T[(^^^^)"'^'c^p|-^(^^-"^^'^)'| • 

V k V k ^ ^^'' ^ 

(58) 



The permutation sum is over all permutations of identical particles (e.g. same spin electrons) and 
the permutation signature e-p = ±1. The initial conditions for Eq. G^ are Wk = 0, m/j = rj^., and 
D = 0. For this ansatz the generator of the norm matrix, Eq. E4, 



N = exp{D + D') J2 ep n[''(^'= + w^k)]"'^' exp {-(nife - m^jV(w^fc + w'vj} 



(59) 



For a periodic system the above equation is also summed over all periodic simulation cell vectors, 
L, with rrifc — m-p^ -^ nife — m-p^ + L. If only the identity permutation is considered the norm 
matrix is easily inverted so that Eq. ^ gives 



8 2 
Wk = ~2wkHD - -w^H^^ 



where 



riife = 


-WkHm^ 


D = 


- 7T" + 

V2 




IdH 


Hqk = 


2dqk ■ 



2y^WiHu 



(60) 
(61) 

(62) 
(63) 



For systems of electrons and ions the full expression for Hq^ and the norm matrix are derived in 
Appendix A. 





FIG. 3. Gaussian approximation for the ground state of a hydrogen molecule for bond length R. The 
top left panel shows the Gaussian mean parameter m for the two electrons. These stay in the center of 
the bond (m = 0) until about R = 2ao and then attach themselves to the separating protons (± R/2). 
The width parameter, displayed in the lower left panel, makes the transition from the optimal value for 
a helium atom, i? = 0, to the hydrogen atom result w = 97r/8an at large R. The right panel shows the 
dissociation energy for the singlet state computed from Eq. p]l (open circles with error bars) and the 
thermodynamic estimator (— dD/d/3) (dashed line) compared to the results of Kolos and Roothan (solid 
line) . 

Application to an isolated hydrogen molecule at low temperature is shown in Figure 0. This 
is for the singlet state (anti-parallel electron spins). The triplet state is considered later after 
a discussion of how to treat permutation terms in the parameter equations. The bond length 
at minimum energy is 1.47 ao, compared with the experimental value of 1.40 ap. The direct 



10 



energy estimator Eq. ^ gives a dissociation energy of 4.50 eV at the minimum compared to the 
experimental value of 4.75 eV. Beyond R — 2, the energy rises quickly toward the value given by 
the Rayleigh-Ritz estimator —dD/d/3. 



VII. ANTISYMMETRY IN THE PARAMETER EQUATIONS 

The determinantal form for the VDM, Eq. pq, is correctly antisymmetric under exchange of 
identical particles. Since ion exchange effects are negligible at the temperatures considered here 
these are ignored. 

The determinantal form leads to A^! terms in the equations of motion for the variational 
parameters presented in appendix A. It was originally hoped that exchange effects could be ignored 
in these equations while retaining the full determinantal form for the VDM but this leads to an 
instability in fermionic systems, e.g. it results in an unphysical strong attraction between two 
hydrogen molecules. 

A practical means of treating all exchange terms, in particular terms involving the potential 
energy, in the variational parameter equations was not found. Instead it was necessary to use 
an approximation similar to that used in the real time computations [|l3| , [16| : only pair exchanges 
in the kinetic energy terms were retained. This will be illustrated for the hydrogen molecule 
after first giving the explicit form for this correction. It is stressed that, unlike the real time 
computations, once the variational parameters are determined the full determinantal form is then 
used in calculating the various averages. 

For two particles with parallel spin, the correction term to the kinetic energy is given by. 



(64) 



AT = — — / dR PAS T PAS - dRpifpj 

PAS ^ Pi{ri)P2{r2) ~ P2{ri)pi{r2) , pi ^ pi{ri)p2{r2) (65) 

Nas = fdR pis , Ni - fdR pj (66) 

For the Gaussian ansatz in Eq. |58| it becomes, 

AT = -^[3(l-w')~Q'] , (67) 



wNq ^ ^ 


— 


w^ 


'WI + W2 , 


w 


=z 



w = wx+W2 , w = „ , = , Q =— (mi-m2) , Nn^w'^e^ -1 . 

2y'WiW2 W 

The corrections to the norm matrix A/" are neglected in order to keep its analytically invertible 
form. The corrections to Hq^ in Eq. p3 are given by 

The correction to dynamics of the parameters follow from Eq. pG to p2l 

Awi = -2 wi (aTd + ^wi AT^, j (70) 

Ami = ~wi AT^, (71) 

AD = -2 {ATd + wi AT^, + W2 AT^, ) . (72) 

These equations lead to an effective repulsion between the Gaussians for two electrons with parallel 
spin if there is significant overlap. As a example of this effect the variational parameters for the 
singlet and triplet states of the hydrogen molecule are compared in Fig. ^. For the triplet state 
parameters the solution including full exchange effects (long dashed line) are compared with those 
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obtained in the kinetic pair exchange approximation (dot-dashed hne). The approximation now 
prevents the Gaussian means for the same spin electrons from collapsing to the bond center at 
lower temperature and is numerically close to the solution for full exchange. 




FIG. 4. Effect of antisymmetry on the density matrix parameters, width and mean, for a hydrogen 
molecule. The protons (large black dots along x axis) are separated by f.Sao and the initial electron 
positions re(/3 = 0) = ±1.5ao along the molecular axis. The solid line for the singlet state (electron spins 
anti-parallel) shows both electrons centered in the molecular bond at low temperatures (large /3). For 
the triplet state (parallel electron spins), long dashed hne the electrons are centered close to the protons. 
The approximation of including only kinetic pair exchanges (dot-dashed line) gives a similar result for the 
mean, with the electrons centered slightly inside the protons but overestimates the Gaussian width (left 
panel) . At high temperature (/3 < 4) exchange is unimportant and the parameters are nearly the same for 
all cases. 

Even at the lowest temperature considered here in the dense hydrogen simulations (5000 K) 
exchange effects between same spin electrons are negligible beyond a few angstroms, i.e. one or 
perhaps two nearest neighbors. Fig. |4| for the triplet state thus overestimates the effect likely in 
dense hydrogen. The main effect of including exchange in the parameter equations is probably to 
prevent the instability mentioned above. 

Fig. a shows an energy comparison for the triplet ground state of the hydrogen molecule. 
First, we compare the Gaussian approximation using only the kinetic exchange term in the 
parameter equations. For the direct estimator, Eq. O, one finds fairly good agreement with the 
quantum chemistry result ||lj]. The thermodynamic estimator gives a somewhat more repulsive 
triplet interaction for R > 2ao. Considering also the Coulomb exchange terms in the Gaussian 
approximation leads to the dot-dashed line for the thermodynamic estimator. We conclude that 
leaving out the Coulomb exchange terms in the parameter equations for efficiency reasons is a 
reasonable approximation in many particle simulations. 
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FIG. 5. Energy of repulsion for the triplet ground state of the hydrogen molecule for bond length 
R. The thermodynamic (dashed line) and the direct estimator, Eq. pll, (circles with error bars) for the 
Gaussian approximation using the kinetic exchange term in the parameter equations are compared with 
the Kolos and Roothan results (solid line) . The thermodynamic estimator for the Gaussian approximation 
with all exchange terms is shown by the dot-dashed line. 



VIII. RESULTS FROM MANY PARTICLE SIMULATIONS 

In this section, we report results from VDM Monte Carlo simulation with 32 pairs of protons 
and electrons in the temperature and density range of 5 000 K < T < 250 000 K and 1.75 < r-g < 4.0. 
Although the Gaussian ansatz VDM will be seen to provide a reasonable model for hydrogen over 
the full density and temperature regime, a large purpose in presenting these results is to serve as 
a base for documenting future improvements from better VDMs and the application of RPIMC. 
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FIG. 6. Proton-proton pair correlation function from VDM (solid line) and RPIMC (dashed lines at 
,=1.75, 2.0, and 4.0 for T < 125 000 K). 



The proton-proton pair correlation functions are shown in Fig. y. For temperatures below 
20 000 K, a peak emerges near IAqq that demonstrates clearly the formation of molecules. The 
comparison with RPIMC simulations |g,|l£l at low density shows that the peak positions agree well 
but RPIMC predicts a significantly bigger height indicating a larger number of molecules. This 
could be explained by the missing correlations in the VDM ansatz. 

At a density of Ts = 2.0, proton-proton pair correlation functions from RPIMC and VDM are 
almost identical. The area under the peak multiplied by the density gives an estimate for the 
molecular fraction. By comparing the estimate for different densities one finds that the molecular 
fraction is diminished when the density is lowered below r^ = 2.0. This effect is well-known and is 
a result of the increased entropy of dissociated molecules. 
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Considerable differences between tlie proton-proton pair correlation functions are found at 
Ts = 1.75 below T = 20 000ii: where VDM shows stiU a fan number of molecules while RPIMC 
predicts a metallic fluid where all bonds are broken as a result of pressure dissociation [p|j2(][| . This 
effect has to be verified by RPIMC simulations with VDM nodes because free particle nodes could 
enhance the transition to a metallic state. 

The peak positions shifts from 1.45ao at a low density of Ts = 4.0 to l.Soo at Vg = 1.75. The 
same trend has been found in the RPIMC simulations [pi but the opposite was reported in i2ll,E3l . 
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FIG. 7. Proton-electron pair correlation functions from VDM (solid line) and RPIMC (dashed lines at 
rs=1.75, 2.0, and 4.0 for T < 125 000 K). 

In the proton-electron pair correlation functions shown in Fig. M, one finds a strong attraction 
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present even at high temperatures such as 250 000 K. At low temperatures, the electrons are bound 
in atoms and molecules. This pair correlation function docs not show a clear distinction between 
the two cases. From studying the height of the peak at the origin multiplied by the density, one 
can estimate the number of bound states at low temperature. Similar to the molecular fraction 
one finds a reduction of bound electrons with decreasing density below r,, — 2.0. The comparison 
with PIMC shows that VDM underestimates the height of the peak. This is probably a result of 
the Gaussian ansatz, which does not satisfy the cusp condition at the proton. 
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FIG. 8. Electron-electron pair correlation function for electron with parallel spin from VDM (solid 
line) and RPIMC (dashed lines at rs = 1.75, 2.0, and 4.0 for T < 125 000 K). 
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FIG. 9. Electron-electron pair correlation function for electron with anti-parallel spin from VDM (solid 
line) and RPIMC (dashed lines at rs=1.75, 2.0, and 4.0 for T < 125 000 K). Note the change in scale in 
the last row. 



Fig. g| shows the effect of the Pauli exclusion principle leading the a strong repulsion for 
electrons in the same spin state. This effect is not present in the interaction of electrons with 
anti-parallel spin (Fig. ^. At high temperature, one observes the effect of the Coulomb repulsion. 
At low temperature, one finds a peak at the origin that is a result of the formation of molecule, 
in which two electrons of opposite spin are localized along the bond. The differences to the PIMC 
graphs can be interpreted as a consequence of different molecular fractions, which has also been 
observed in Fig. |6[ 
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FIG. 10. Average width of the Gaussian single particle density matrices as a function of temperature 
for different densities 

The average width w of the Gaussian is shown in Fig. ^ as a function temperature and 
density. At high temperature and low density, one finds only small deviations from the free particle 
limit. These become more significant with increasing density and decreasing temperature. At low 
temperature, the attraction to the protons dominates, which leads to a decreasing average width. 
Finally bound states form and the width approaches a finite limit. At low densities, this is close 
to the ground state width of the isolated molecule 3.138 ag. 
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In Fig. 



the direct estimator 



11 , we compare the internal energy from the thermodynamic estimator in Eq. |5^ and 
Both agree fairly well at low density. Differences build up with increasing 
density and decreasing temperature. Comparing with RPIMC simulations, one finds that the VDM 
energies are generally too high. The magnitude of this discrepancy shows the same dependence on 
density and temperature like the difference between the two VDM estimators. The difference to 
the RPIMC results could be explained by the missing correlation effects in the VDM method. 

At high temperature, the thermodynamic estimator always gives lower energies than the direct 
estimator. Below T — 2h 000 K, the ordering is reversed. This is consistent with the results from 
the isolated atom and molecule. The consequence is that the direct estimator is actually closer to 
the value expected from RPIMC simulations. How ever, it should be noted that this estimator is 
not thermodynamically consistent (see section V B ) . 
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FIG. 12. Pressure versus temperature in high and low temperature range. VDM pressure is calculated 
from virial relation using both the direct and thermodynamic estimators for kinetic and potential energy. 

In Fig. n2l we compare pressure as a function of temperature and density from the two 
VDM estimators with RPIMC resuhs. At low density, the agreement is remarkably good. With 
increasing density and decreasing temperature, the difference grows. For densities over r^ = 2.0 
below 10 000 K, one finds a significant drop in the direct estimator for the pressure. We interpret 
this effect as a result of the thermodynamic inconsistency. 
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FIG. 13. Comparison of experimental and several theoretical Hugoniot functions 



Fig. 03, compares the Hugoniot from Laser shock wave experiments p3,p4 with results 
from several theoretical approaches (Sesame data base by Kerley |2^] (thin solid line), linear 
mixing model by Ross (dashed line) [E6[, tight- binding molecular dynamics by Lenosky et.al. 
pTI (dash-dotted line), Fade approximation in the chemical picture by Ebeling et.al. ^^ (dotted 
line), RPIMC simulations |p^ (triangles), VDM direct estimator (full diamonds) and VDM 
thermodynamic estimator (full circles)). The long dashed line indicates the theoretical high 
pressure limit p = Apo of the fully dissociated non-interacting plasma. In the experiments, a 
shock wave propagates through a sample of precompressed liquid deuterium characterized by its 
initial state, {Eq, Vq, po). Assuming an ideal shock front, the variables of the shocked material 
{E, V, p) satisfy the Hugoniot relation Q, 



1 



H = E-Eo + -{V-Vo){p + Po)^0 



(73) 



The initial conditions in the experiment were T = 19.6 K and p — 0.171 g/cui^. We set Vq = 39.1 A'^ 
and Po ~ 0. We show two VDM curves based on the thermodynamic and direct estimators. 
For £^0, we use the corresponding value of the ground state of the isolated hydrogen molecule, 
£*'» = -0.955 Ha and E^'"' = -1.124 Ha. 

We expect the difference of the two estimators to give a rough estimate of the accuracy of the 
VDM approach. At high temperature, the difference is relatively small and agreement with RPIMC 
simulations is reasonable. Both VDM estimators indicate that there is maximal compressibility 
around 1.5 Mbar. However, in this regime of high density and relatively low temperature a more 
careful study seems unavoidable. We suggest RPIMC simulations using the VDM nodal surface to 
restrict the paths. 
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IX. CONCLUSIONS 

The VDM approach provides a way to systematicaUy improve the many particle density matrix. 
Already the simplest ansatz using one Gaussian to describe the single particle density matrices gives 
a good description of hydrogen in the discussed range of temperature and density. The method 
includes the correct high temperature behavior and shows the expected formation of atoms and 
molecules. The thermodynamic variables are in reasonable agreement with RPIMC simulations. 
The presented Gaussian ansatz can be improved in several ways. One could use a sum of Gaussians, 
add underestimated correlation effects by including a Jastrow factor in the ansatz or use a two-step 
path integral. Further one can use this essentially analytic density matrix to furnish the nodal 
surface in RPIMC simulations, replacing the free particle nodes by a density matrix that already 
includes the principle physical effects. This level of accuracy seems to be required to determine a 
Hugoniot function that is very sensitive to the different level of approximations made by various 
theories. 
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APPENDIX A: GAUSSIAN APPROXIMATION INTERACTION TERMS 

The general equations for the variational parameters g in a parameterized density matrix, from 



Eq. 21 



are 

IdH ^ ■ 

where 

H= I pHp dR = / pHpidR (A2) 

and the norm matrix 

Nji = / Pj p^ p^ dR = lim ——J (A3) 



with 



N = J p(R, q; (3) p(R, q' ; /3) dR . (A4) 



The subscript / in Eq. A2 indicates that only one p needs to be antisymmetric and the identity 
permutation can be used in the other. (We are also dropping 1/A^! prefactors which are the same 
for the norm matrix and thus cancel out.) This appendix contains the detailed formulae for these 
equations for a parameterized Gaussian density matrix applied to a Coulomb system. 

Repeating Eq. |5^ the parameterized variational density matrix is an anti-symmetrized product 
of one-particle density matrices, 

p{B.,R',(3)^Y.'r\{pi{vkyv,^P)^Y.'ve''\{{^w-pX'''^e^v{~ — {vk-rav,f\ (A5) 
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where the amphtude D and the widths Wk and means nife are the variational parameters. The 
permutation sum is over all permutations of identical particles (e.g. same spin electrons) and 
e-p = ±lis the permutation signature. The initial conditions arc Wk = 0, irife = r'^,, and D = Q. 
For this ansatz the generator of the norm matrix, 

N ^Y^ev [][7rK + ^pjr'^' exp {-(mfe - m^jVK + <J} exp(D + D') . (A6) 

V k 

For a periodic system the above equation also is summed over all periodic simulation cell vectors, 
L, with rrifc — irip -^ m.k — m'^ + L. Using this the components of the norm matrix are then: 



Mdd = ^ e-pN-p 



A/m.D = y^e-p 



^f^ 



WiD 



E 



e-p 



— 2(mi — m-p.) 


Wi + W-p. 

-1 ^ 


■3 


Wi + W-p^ ) 


2 



Np 



3 ( m, - mpj^ 

Wj + Wp^ 



Nv 



A/'miirij = 2^ ep 



2^j-,P. y ^ ^ (m,~mp,)("^J-"^p-0 

Wi + Wj {wi + WpJ {Wj + Wp-i) 



TVp 






V 



Oj.V, 



(mj - m^ 



Wi + Wj {wj + Wj,-i) \ 2 {wj + Wj,-i) 



2 (mi -mpj 

Wi + Wp. 



(A7) 
(A8) 
(A9) 

(AlO) 

Np, (All) 



E 



ep 



^i,^. 



(Wi + Wp, 



3 2(mj-mpj2 
2 



Wj + Wp^. 



1 



3 _ (mi - mpj^ 
2 



Wi + WPj 



3 (mj-mp-i)2i 



Wj + Wp-l 



(Wi + WpJ{Wj + Wp-l) 

Np 



where 



A^P = e^^ n '''''' 1 ("'^■+"'-.) 
J (7r(w, +wpJ)3/2 



A^i 



p-i 



The Hamiltonian for a periodic system of electrons and ions 

w = -2E ^?+EE ^(^^^■) - E E ^/^(r^^) + E ^^t-'i + ^" 

2—1 i<j i I i 

where the purely ionic terms 

KI' I 



(A12) 



(A13) 



(A14) 



(A15) 



The Ewald potential, V'(r), which includes interactions with periodic images and incorporates 
charge neutrality. 



*W ^ E =^^1^1,^ + E tS =P<-^V4GV ^ ^ S i&=p(.k . 



k^^iO 



f7fc2 






(A16) 



where f2 is the periodic cell volume and G an arbitrary constant. The Madelung term in 7i is 
the interaction energy of an electron with it's periodic images and neutralizing background (e.g. 
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UMad = — 1.41865/L for a simple cubic simulation cell, the usual case). To do the integrals we 
represent the Gaussians by their Fourier series 



,|_,3,,^.-.,_-.,^^l, 



'fc^-!i)/8 ik-(r-m) 



(A17) 



L k 

and in the interaction terms use the Fourier representation for ^j{v). This finally gives 

H = ^€v{Kv + Uv}Nv (A18) 



with 



3 (mi-m-pi)2 



Wi+W-pi {Wi+W-piY 



(A19) 



i<j i 1 i 

where Wi = WiW-pi/ [wi + w-pi) and rhi = (niiW-pi + m-piWi) / {wi + tw-pi) . The interaction integral 

47r _fc^ 



iy(r, w) = Y^ ^e-'^e^-'' 



(A21) 



fe^^iO 



VF is symmetric in r when the periodic cell has inversion symmetry. Continuing, the left hand side 
of Eq. ^ is 

IdH 



^^ = 2aD-^ 



IdH Iv- { ,dKv dUv.^r 



H„ 



1 dH 1 
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2 am, 2 
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(A22) 
(A23) 

(A24) 

(A25) 
(A26) 
(A27) 
(A28) 



where we have used the fact that terms in Vi and V ^i give the same contribution under the 
permutation sum and so combined them. The derivatives of the interaction integral are. 



dUp 2wpi 



dnii Wi + wpi 
dUp '2w-pi 



dwi {wi + wviY 



Y^ will (m, - m^,w, + Wj)-Y^ Z/W[il {m, - R/, ', 



^■p^ Yl ^'^' ("^' " "^j' *' +^j)-J2 ^^^'^' ("^' " ^^' *') 



(A29) 



,j#' 



+ (mp,-mO • I ^wW(m, -mj,w, +%)-^Z/WW(m, -R/, 



(A30) 
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where WWand W^'^^ denot e the derivatives of W with the first and second argument. Comparing 
equation A21 and Eq. Af 6 the interaction integral may be written as 



Wir,w)^^ir)-J2 



erfc 



k+L| 



|r + L| 



TTW 



and its derivatives as: 






r + L 

r + LP 
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/TTW 



exp(-|r + L|Vw) 



^ exp(-|r + L|Vw) n 



■3/2, 



Vl 



For an isolated system (L — > oo) and these would simplify to, 

erf \rl^/w ] 



will (r, w) = - ^ ( erf [r/ V^ ] - 



2r 



,-r^h- 



W^'^\v,w) 



1 



-r I w 



Wy/TTW 

At /3 = the initial derivatives for the variational parameters reduce to 

Wi = 2 



(A31) 

(A32) 
(A33) 

(A34) 
(A35) 
(A36) 



(A37) 
(A38) 
(A39) 



For large numbers of electrons it is not possible to treat all permutations. Here the 
approximation discussed in section VII is used where the kinetic pair exchange corrections given 
there are added to the identity permutation term derived here. 
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